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Abstract 


The hypersonic viscous flow over a flat plate at zero angle of attack and wedge is de- 
termined by solving Navier-Stokes equations, using a time-marching technique. The 
steady state solution is attained asymptotically. The computation is performed in 
a rectangular domain. The method used for solving these equations is the explicit 
MacCormack finite difference scheme . This is a predictor-corrector scheme which is 
second order accurate. Shock-capturing method has been used where the location of 
the shock is obtained as a steady state solution without having an apriori knowledge 
about it. A code has been developed which generates numerical solution for any given 
Mach number. Detailed results for Mach numbers 4, 8 and 12 are presented and com- 
pared with existing results. The results obtained in this case have also been compared 
with those obtained from boundary layer equations and Newtonion approximations 
which shows that the coefficient of pressure on the body comes closer to the Newto- 
nian approximation with increase in Mach number . Effect of number of grid points in 
domain on the result has been analysed. Results from both uniform and non-uniform 
grids have been obtained and compared. Reduction in coefficient of pressure and 
friction coefficient with Mach number has been observed. In addition to this it has 
been observed that the rate of increase in shock-wave strength for a wedge is higher 
than that for a flat plate. It has been observed from the results obtained that the 
strength of the shock wave increases and shock layer thickness decreases with increase 
in Mach number, as expected. In addition to this increase in strength of shock wave 
with increase in wedge angle has also been observed. 
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Chapter 1 


Introduction 


1.1 Special Features of Hypersonic Flows 


According to modern classification flow regimes are divided into sub-sonic flows (upto 
Mach Number 0.8), transonic flows (from Mach Number 0.8 to 1.2), supersonic flows 
(from Mach Number 1.2 to 5) and hypersonic flows (above Mach Number 5). Lin- 
earized theory can be used for studying subsonic and supersonic flows but the study 
of hypersonic flows is complicated, and requires special treatment because of the fol- 
lowing reasons: 


1. Shock layer is very thin for hypersonic flows. This thin shock layer can create 
some physical complications, such as the merging of the shock layer itself with 
a thick, viscous boundary layer growing from body surface. The shock-wave 
boundary-layer interaction becomes important at low Reynolds numbers. 

2. Strong entropy gradients are generated in the nose region and the boundary 
layer along the surface grows inside this entropy layer and is affected by it. This 
entropy layer causes analytical problem due to simultaneous effect of boundary 
layer and entropy layer in the vicinity of the wall, while performing a standard 
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boundary layer calculation on the surface. 

3. Boundary layer thickness in hypersonic flow is higher which can exert a major 
displacement effect on the inviscid flow outside the boundary layer. Due to the 
extreme thickness of the boundary layer flow the outer inviscid flow is greatly 
changed which in turn affects the growth of the boundary layer. This major 
interaction between the boundary layer and the outer inviscid flow is called 
viscous interaction. Viscous interaction affects the surface pressure distribution, 
hence lift, drag and stability of hypersonic vehicles. Moreover, skin friction and 
heat transfer are increased by viscous interaction. 

The boundary layer on a hypersonic vehicle is so thick that it merges with shock 
■wave. In this case the shock layer is treated as fully viscous and the conventional 
boundary layer analysis is abandoned. 

4. Boundary layer and normal and nearly normal shock in the nose region cause 
high temperature in hypersonic flow, which makes the flow chemically reacting. 
High temperature chemically reacting flow's can influence lift, drag, and moments 
on a hypersonic vehicle. However the most important dominant aspect of high 
temperature is the resultant high heat transfer rates to the surface. Aerodynamic 
heating dominates the design of the hypersonic vehicle. 

In addition, high-temperature causes ionization producing free electrons w r hich 
absorb radio frequency radiation. This results in communication blackout i.e. 
no transmission of radio waves either to or from the vehicle. Therefore, the 
accurate prediction of electron density within the flow' field is important. 

5. Since some hypersonic vehicles, fly at or through the outer region of atmosphere 
they experience the effect of low density. With increase in altitude mean free 
path (A) and similarity parameter governing density regime (Knudsen number 
Kn = A/L), increase. Continuum Navier-Stokes equations are applied up to 
Kn<0.2. Slip effects are included in these equations when Kn>0.03. The effect 
of free molecular flow begins around a value of Kn=l and extends out to the 
value of Kn=oo. Hence the transitional regime is essentially contained between 
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0.03<Kn<l. 


1.2 Solution Methods of Hypersonic Viscous Flow 
Problems 


Various approaches have been applied for the solution of hypersonic viscous flow prob- 
lem. Before the advent of computational fluid dynamics exact solutions of complete 
Navier-Stokes equations for practical problems were difficult. This is the reason of 
simpler viscous flow solutions. The approximations used at that time involved sepa- 
rate calculations of the boundary layer and outer inviscid flow, and then a coupling 
of these separate calculations to take into account the viscous interaction. The weak 
and strong interaction serve a useful purpose in providing convenient correlations and 
predictions based on approximate theory. By a suitable order of magnitude reduction 
of Navier-Stokes equations boundary layer equations are obtained. The compressible 
boundary layer equations are the same for all flow regimes i.e subsonic, supersonic 
and hypersonic (neglecting high temperature effects). Various methods have been 
introduced to solve boundary layer equations. 


1. The difference- differential method gives an exact solution of the general boundary 
layer equation. Smith and Clutter [1] have successfully solved these equations. 
They have stated that meeting boundary conditions efficiently has been the most 
difficult part of the entire flow problem. 

2. Finite difference method was suggested by Blottner [2], these are inherently faster 
and more accurate than the other methods. He solved these equations for multi 
component flow with finite chemical reaction. DiCristina has obtained experi- 
mental results for a sharp 8 degree cone at Mach Number 10 [3] . Stetson has 
shown the results at Mach Number 6 over a cone [4]. 
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If we consider the entire flow as fully viscous and no partition between boundary 
and inviscid flow 7 , the viscous interaction effect can be calculated exactly. 

There are basically three approaches used to solve the fully viscous flow problem 
which provide every information about the flow such as shock shape, detailed flow 7 
variables between shock and body, skin friction, heat transfer, lift, drag, moments etc. 


1. Viscous Shock Layer Technique 

This technique was introduced by Davis [5] who derived a set of equations w'hich 
approximate the Navier-Stokes equations and solved them for a viscous shock 
layer over a blunt body at hypersonic speed. An order of magnitude reduc- 
tion was applied to the Navier-Stokes equations. The terms are kept in second 
order in 1/y/Re (Re = Reynolds number) to obtain the governing equations. 
Since these equations are parabolic hence downstream marching finite difference 
method is used. This method has found wide application to chemically reacting 
viscous flow. 

2. Parabolized Navier-Stokes Solutions 

This approach uses equations which contain more terms than the the viscous 
shock layer equations and hence more accurate but are not as accurate as the full 
Navier-Stokes equations. They are obtained from full Navier-Stokes equations 
(w T ith all time derivatives zero) simply by neglecting all viscous terms w r hich 
involve derivatives in the streamwise direction. 

With two exceptions, these equations are a mixed syst em of parabolic-hyperbolic 
partial differential equations and hence are solved by a downstream marching 
procedure. The two exceptions are, 

(a) In subsonic region, the pressure is constant in the direction normal to, 
the surface, equal to its value at the first grid point at which supersonic 
flow 7 exists. This is done to preserve the parabolic nature of Parabolized 
Navier-Stokes (PNS) equations. 
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(b) The volumetric heating term, namely, can destroy the parabolic behaviour 
of the system. For flows where such volumetric heating does not occur, 
such a problem does not exist. 

McWherler [6] has solved parabolized Navier-Stokes equations successfully in his 
work. He has calculated the flow over slender blunt nose cone at small angle of 
attack and compared the results with boundary layer solutions. 

3. Full Navier-Stokes Solution 

In this approach all the terms present in the Navier-Stokes equations are con- 
sidered and is the ultimate viscous flow calculation. Hence this is used to obtain 
solution over a flat plate and a wedge. Navier-Stokes equations are a system 
of partial differential equations with a somewhat mixed hyperbolic, parabolic 
and elliptic behaviour. Due to their elliptic nature these equations cannot be 
solved by a downstream marching technique. Hence full Navier-Stokes equa- 
tions have been solved by using a time dependent finite difference method. The 
major reason for using the time dependent method is that resulting unsteady 
Navier-Stokes equations are a mixed set of hyperbolic-parabolic equations for 
all subsonic and hypersonic flows. As a result, a very complicated flow field can 
be calculated as an initial value problem. An additional advantage is that since 
the complete Navier-Stokes equations are solved the shock wave, shear layer and 
a wall boundary layer are captured automatically in the solution without prior 
knowledge of their location or even existence. 


1.3 Methods Available for Solution of Full Navier- 
Stokes Equations 


Many numerical approaches are available for the solution of full Navier-Stokes equa- 
tions using both implicit and explicit finite difference methods. Nearly all the methods 
are second order accurate in space and are first or second order accurate in time. The 
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explicit finite difference methods available are as follows: 


1. MacCormack Method:- This method was introduced by MacCormack [7]. It uses 
forward difference for all special derivatives in predictor steps and backward 
difference in corrector steps. This eliminates any error due to one sided differ- 
encing. The method and it’s application to the Navier-Stokes equations is given 
in Chapter 2. 

2. Hopscotch Method'At is an unconditionally stable method. It involves two sweeps 
through the mesh. Let i and j denote the number of grid points along x and y 
axis respectively and n be the time-step under consideration. For the first sweep 
the variables in the next time step e.g.: Uij 1 is computed at each grid point ( for 
which i+j+n ) is given by the simple explicit scheme and for the second sweep 
these variables are computed at each grid point( for which i+j+n ) is odd by the 
simple implicit scheme. 

3. Leap Frog/Dufort -Frankel Method 

This is a one step scheme and is first order accurate. In this method, time 
derivative is given by central-difference as shown below 

(duY = tff 1 - 

\dt) . 2 At 

First order special derivative is given by central difference. This scheme is better 
suited for calculation of steady solutions, where time accuracy is not important. 

4. Allen-Cheng Method 

This scheme is a first order accurate, predictor-corrector scheme. In predictor 
step, time derivative is given by the difference of predicted value in the next time 
step and value at the present time, while in corrector step it involves difference 
of corrected value in the next time step and the present value. This method 
allows larger time steps when viscosity is large. 
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5. Lax-Wendroff Method 


It is a second order accurate two step method. In the first step variables at the 
half of the next time are determined. In the second step the variables obtained 
in the first step are used to find them at the next time step. 

When the methods other than MacCormack method are applied to the compli- 
cated compressible Navier- Stokes equations certain problems can arise. For example 
the mixed derivative terms can create a problem for the Hopscotch method. If these 
terms are differenced in the usual manner the Hopscotch method is no longer explicit 
since a matrix inversion is required. All the above methods except the Lax-Wendroff 
scheme are first order accurate in time and so they cannot be used to accurately 
calculate the time evaluation of the time field. Being one of the better methods Mac- 
Cormack method is used to solve the problem. 


1.4 Assessment of Existing Literature 


Initially approximate solutions were derived for hypersonic flows .For a tangent cone, 
one pressure approximation was given by Hays and Probstein [8] using constant density 
approximation and Maurice [9] derived a simple, approximate solution for a cone. Datv 
and Maurice [10] obtained an approximate analytical solution by means of constant 
density approximation, for hypersonic flow past an inclined cylinder. 

As stated earlier, there are various approaches which constitute the hypersonic 
viscous flow problem, but the solution of full Navier-Stokes equation is the exact 
solution. An attempt has been made in the succeeding paragraphs to identify the 
major steps and advances in the field of hypersonic viscous flow problem. 

Work in this field started with the introduction of various methods to solve the 
boundary layer equation, since before the advent of computational fluid dynamics it 
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was not possible to solve the full Navier-Stokes equations. The pioneering work of 
Moretti and Abett[ll] gave a time dependent computational method for blunt-body 
flows. The computation was performed for shock layer only, assuming shock layer 
is a discontinuous surface, dependent on time. Results for 2-D blunt bodies and for 
axisvmmetrie ones were obtained and compared with experimental data. 

As an extension of his work, MacCormack [12] gave a numerical method for solving 
the equations of compressible viscous flows. This method is second order accurate in 
space and time, unconditionally stable, preserves conservation form, requires no block 
or scalar tridiagonal inversions, and thus more efficient than the existing methods. 
The validity of the method was checked by externally generating shock wave incident 
upon boundary layer on a flat plate. 

After introduction of viscous shock layer technique by Davis[5] Horstman [13] 
gave solution for hypersonic/turbulent boundary layer interaction flows. Lind and 
Lewis [14] studied unsteady characteristics of shock-wave interaction (for which the 
peak pressure, heat transfer rate and surface pressure distributions are sensitive to 
upstream thermodynamic flow conditions, shock length and Mach number) at Mach 
Number 8 using a time accurate scheme for solving the thin layer approximation to 
the two dimensional Navier-Stokes equations. They showed by calculations that the 
peak surface pressure, the impingement location of the supersonic jet and the time 
required for the development of the interaction are strong functions of the impinging 
shock, location. For a range of impinging shock location of interaction is found to be 
unsteady. 

Gupta et al. [15] presented numerical solutions from the steady viscous shock-layer 
equations for the hypersonic laminar and turbulent flow of a perfect gas over long slen- 
der bodies. Kussoy and Horstman[16] presented experimental data for two three di- 
mensional intersecting shock wave/turbulent boundary layer interaction flow’s at Mach 
number 8.3. Gupta et al. [17] presented a solution procedure that improved the com- 
putational efficiency of the viscous shock layer technique, especially for long slender 
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bodies. Cheatwood and DeJarnette [18] had developed an axisymmetric approximate 
method which could calculate laminar (perfect gas, non-equilibrium, equilibrium) and 
turbulent(perfect gas and equilibrium) viscous hypersonic flow over blunt nosed body. 
Simplified form of viscous shock layer equations were used. Results were obtained for 
flow over analytic body shapes. Surface heat transfer and pressure predictions were 
generally in good agreement with viscous shock layer results. 

Bigdeli et.al. [19] used turbulent viscous interaction theory to understand the hy- 
personic flow past an expansion corner. The paper of Luca and Cardone [20] dealt with 
an experimental investigation carried out to study some aspects of shock/boundary 
layer interactions in normal two dimensional hypersonic wedge flow' over flat-plate/ramp 
configuration. 

The parabolized Navier-Stokes (PNS) equations have also been used to successfully 
compute the 3-D hypersonic viscous flow over variety of body shapes at angle of attack. 
In this series Tannehil et.al. [21] had developed a parabolized Navier-Stokes code to 
compute the steady supersonic viscous flow around arbitrary body shot at high angle 
of attack. A nonorthogonal 3-D co-ordinate frame was employed which permitted the 
code to march with solution surfaces. The code w r as used to calculate the laminar flow r 
over a slab delta w r ing with 70 degree sweep at angles of attacks upto 41.5 degrees 
and Mach Numbers of 6.8 and 9.6. Barnett and Davis [22] described a procedure 
for the calculation of strong viscous inviscid interactions in 2-D laminar supersonic 
interactive flow’s with and without separation. The equations solved were parabolized 
Navier-Stokes equations. Results were obtained for flow’ past flat plate related bodies 


Dash et.al. [23] described parabolized Navier-Stokes models under development 
for the analysis of 3-D supersonic and subsonic propulsive jet mixing problem, and 
Kimmel et.al. [24] assumed the stability of boundary layers on sharp nosed cones 
with elliptical cross sections using linear stability theory and cross flow correlations. 
PNS computer codes were used to calculate the mean flow about cones with various 
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eccentricities at a free stream Mach Number 7.95 and a Reynolds Number 3.3Xl0 6 ra 1 . 

For complete and exact solution of hypersonic viscous flow problem full Navier- 
Stokes equations have been solved. In this context Kumar and Salas [25] obtained 
Euler and Navier-Stokes solutions of the supersonic shear flow past a circular cylinder. 
The results obtained indicate that the Euler equations may correctly predict certain 
high Reynolds number separation phenomena in flows with a natural inviscid vor- 
ticity source. Hung and Timothy numerically solved the compressible Navier-Stokes 
equations for hypersonic flow over a 3-D ramp with a narrow expansion slot. They 
had shown that wall cooling reduces the thickness of boundary layer and hence in- 
creases the flow expansion substantially. Hung and Barth [26] solved compressible 
Navier-Stokes equations for hypersonic flow over a three dimensional ramp with a 
narrow expansion slot. Solutions were obtained for various wall temperatures and slot 
widths. 

For the very first time Shang and Scherr[27] presented the numerical simulation 
of the flow-field around a complete lifting body configuration utilizing the Reynolds 
averaged Navier-Stokes equations. The numerical solutions generated for the experi- 
mental aircraft X24C-10D, at a Mach Number of 5.95, exhibited accurate prediction 
of detailed flow properties and integrated aerodynamic coefficients els well. Zoby et.al. 
[28] in their technical notes gave an idea about hypersonic non-equilibrium viscous 
solutions over slender bodies. 

Oberkampf and Aeschliman [29] presented aerodynamic force and moment mea- 
surements and flow visualization results for a hypersonic vehicle configuration at Mach 
Numbre 8 and after that Walker[30] developed parabolized and Navier-Stokes code to 
predict flow field solution around a hypersonic vehicle. Aerodynamic forces and mo- 
ment predictions from the codes are compared with experimental data obtained in 
Ref. [29] . Lockman et.al. [31] established a benchmark experimental database for a 
generic hypersonic vehicle shape for validation and /or calibration of advanced compu- 
tational fluid dynamics computer codes. They included results from the comprehensive 
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test program conducted in the NASA Ames 3.5 ft wind tunnel for a generic all-body 
hypersonic aircraft model. Computational results obtained were compared with exper- 
imental results. Papuceuoglu [32] experimentally investigated the high heat transfer 
rates and flow patterns in the viscinity of 90 degree axial corner at a Mach Number 
of 6 and a Reynolds number range from 7 to 22.5X10 6 . 

Tai and Kao [33] developed a full Navier-Stokes code for predicting the aerody- 
namic properties of slender bodies. An explicit upwind flux-difference split scheme 
combined with a multistage method was implemented for solving the axisymmetric 
full Navier-Stokes equation. Experimental data were found to be in good agreement 
with the code predictions. 

Apart from the works stated above some other works have also been done in the 
field of hypersonic viscous flow, e.g.: Kim et.al. [34] obtained numerical solutions for 
the Navier-Stokes equations for laminar hypersonic flow' around a flat based projectile 
using perfect gas, equilibrium air and non-equilibrium air models. The effects of 
different gas models were investigated for a free stream condition of Mach Number 
20 and at an altitude of 73 km. Both implicit and explicit integration schemes were 
used. Morgenstern and Chokami [35] solved time dependent compressible Navier- 
Stokes equations for hypersonic flow over a cavity. 

Mallett et.al. [36] implemented a kinetic theory based Navier-Stokes solver on a 
parallel super computer to study the leeward flow-field of a blunt nosed delta wing at 
30 degree incidence at hypersonic speed. Computational results were presented for a 
series of grids for both inviscid and laminar viscous flow. Kopriva [37] presented steady 
solutions of high speed viscous flow over blunt bodies. The region within a shock 
layer was divided into subdomains so that internal layers can be well resolved. At 
interface between subdomains, the advective terms were upwinded and viscous terms 
were treated by a penalty method. The method was applied to various hypersonic 
flows. Results were compared with the experimental data and to a finite difference 
solution. 
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Habashi et.al. [38] presented a finite element, segregated method for 2-D hy- 
personic thermochemical non-equilibrium flow, with emphasis on efficiently resolving 
shock- wave by mesh adaptation. The governing equations were decoupled into three 
systems of partial differential equations, gas dynamic, chemical and vibrational sys- 
tems and then solved in a sequential manner. 


1.5 Scope of the Work 


The results obtained in the work give the exact value of skin-friction coefficient and 
the rate of heat transfer at the surface, which are useful parameters for the selection 
of materials, design, calculation of thrust and rate of cooling required for hypersonic 
vehicles (e.g. space vehicles, missiles etc.). 

Hence the work done is very useful in the space and missile programme. Results 
may be helpful in the exact calculation of the bluntness required for the nose of a 
missile or spacecraft (which decides the temperature gradient in that region) and 
optimize the drag. The problem finds its important application in the analysis of 
effect of aerothermodynamics of the flow on the structural dynamics of the body. As 
observed at higher Mach numbers temperature and pressure on the body are very 
high. Since high temperature reduces the stiffness of the structure, both temperature 
and pressure at hypersonic speed causes large deflection on the body (wings of the 
aircraft). Thus limit of the speed can exactly be defined by considering this effect. 

The present problem can be extended to a three-dimensional problem by consid- 
ering the third spatial parameter of the Navier-Stokes equation. 
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Chapter 2 


Problem Statement and Governing 
Equations 


2.1 Problem Statement 


Consider a viscous fluid flowing steadily at a hypersonic speed over a wedge of angle 
a. To predict the macroscopic behavior a domain has been defined. In flow boundary 
of the domain is at leading edge of the wedge. The lower boundary is at the surface 
of the wedge while upper boundary is at a distance of 55 (5=boundary layer thickness 
and is given by the expression 5 = ^==, where L is length, of plate) from the wedge 
and parallel to it. The problem is confined to the face of the wedge and does not 
include backward facing step at the trailing edge of the wedge. 

A cartesian coordinate system is chosen, x and y distances are measured by 
considering leading edge of the wedge as origin, wedge as x axis and inflow boundary 
as y axis. Flow approaches with a uniform velocity V, x The components of velocity 
in x and y directions are u and v respectively. It is assumed that the flow is two 
dimensional. 
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2.2 Governing Equations 


2.2.1 Continuity Equation 


The conservation of mass applied to a prefect gas passing through an infinitesimal, 
fixed control volume yields the following equation of continuity 

f + V.(pV) = 0 (2.1) 


For a cartesian coordinate system, where u, v, w represent the x, y, z components, 
respectively, of velocity vector V and p represents density, the above equation becomes 


dp , d(pu) djpv) djpw) 

dt ^ dx dy dz 


( 2 . 2 ) 


2.2.2 Momentum Equation 

Newton's second Saw applied to a prefect gas passing through an infinitesimal, fixed 
control volume yields the following momentum equation 

^ + vpW = pf + v-n ij (2.3) 

where pf is the body force per unit volume which is mainly the gravitational force, 
hence 

pf = pg (2.4) 

and v • Hy represents the surface forces per unit volume. In this equation pW is a 
tensor so that y • pVV is not a simple divergence.This can be explained as 

V • pVV = pV • yV + V (v • pV) (2.5) 

Putting expression for y • pVV in Equation 2.5 and simplification by using the con- 
tinuity equation gives 

DV 

P ~Dt = pf + V ' ( 2 -6) 
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where pf is the body force per unit volume which is mainly the gravitational force. 


These forces are applied by external stresses on the fluid element. The stresses 
consist of normal stresses and shearing stresses and are represented by the components 
of the stress tensor Iljj, which can be expressed in compact tensor notation as 

n >j = ~Pkj + P + Si ^'faT k = 2 ’ 3) 

where p is pressure and Si j is Kronecker delta function ( 6 t j = 1 if i=j and Sij = 0 if i^ 
j);ui, u 2 , U 3 represent the three components of the velocity vector l r ;rc 1 , rr 2 , x 3 represent 
the three components of the position vector;// is the coefficient of dynamic viscosity 
and p' is the second coefficient of viscosity.// and p! are related by the expression 

K=\, u-tV (2.8) 


where K is the coefficient of bulk viscosity and is negligibly small except in the study 
of the structure of shock waves and in absorption and attenuation of acoustic waves, 
hence 

p' = ~\p (2-9) 

Thus equation can be written as 


DV 8 

f -DT = ft -™+gT i 


( dui du t \ 2 . duk 

^ + dxi j 3 dl ’ j ^dx k 


( 2 . 10 ) 


or 

^ = pf-VP + J-M (2.11) 

where is the viscous stress tensor. 


Writing equation for a cartesian coordinate system, obtaining them in conservation 
law form and assigning Ui = u, u 2 = v, u 3 = w gives 


d(pu) + d_ 

dx 


dt 


d(pv) d 


+ ^ (p « 2 + P~ t x . x ) + ~{puv~ r xy ) + — (puw - r xz ) = pf x 
0/o \ d 


+ — ( PVV - T xy ) + — (pu 2 + P - Tyy) + — ( PUW - T yz ) = pfy (2.12) 


dt ' dx ' r "~ ‘ dy 

d(pw) d 0 . 0/ 2 \ , 

+ (PUW - T xz ) + — (pVW -T yz ) + — (pl/T + p - T ZZ J = pf x 


dt 
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where the components of the viscous stress tensor Tij are given by 


2 ( du dv dw\ 

Txx 3 ^ ^ Q x Qy Q z J 

2 f n dv du dw\ 

Tyy = 

2 ( n dw du dv\ 

Tzz 3^ \ dz dx dyj 

( du dv \ 

Txy ~ tl \dy + !h)~ Tyx 
( dw du \ 

T ~ = " + sy) = T " 

( dv dw\ 

T »* = '‘l& + 8Fj = r " 

(2.13) 


f = / I i + / i ,j + /.k (2.14) 

where i, j and k are unit vectors along x, v and z directions respectively. 


2.2.3 Energy Equation 

The first law of thermodynamics applied to a prefect gas passing through an infinites- 
imal, fixed control volume yields the following energy equation: 

^ + v£.v = ^- vq+pf-v+v (n,j -V) (2.15) 

where E t is the total energy per unit volume given by 

E t = p + PotentialEnergy + • • (2-16) 

and e is the internal energy per unit mass. ^ in Equation (2.15) represents the 
rate of increase of the total energy per unit volume in the control volume while V-Q 
represents the rate of total energy lost by convection (per unit volume) through the 
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control surface. For cartesian coordinate system equation becomes 

~ “ It ~ p ( fxU + fyV + fzW ) + + pu ~ UTxx ~ VTxy _ ™ Tr2 + 9l) 

+ ^ (£)*-’ + pu - ur Iy - ur yy - wT yz + q y ) + {E t w + pw — ur xz — vr yz — wt zz + q z ) — 0 
^ 2 ( 2 - 17 ) 

where 

q = g T i + g y j + g*k. ( 2 - 18 ) 

In terms of heat conduction 



2.3 Vector Form of Equations 


For consistent application of finite difference algorithm to the governing fluid dynamic 
equations, they are combined in compact form. The compressible Navier-Stokes equa- 
tions in two dimensional cartesian coordinates without body force or external heat 
addition is expressed as 


dU ^F 

dt dx + dy 


( 2 . 19 ) 


where U, E and F are vectors given by 


P 


U = 


pu 

pv 


E t 
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pu 

pu 1 2 + p — T xx 
pUV - T xy 

(E t +p)u- UT XX - vr xy 4- q x _ 

pv 

PUV + P~T X y 
PV 2 +P~Tyy 

(E t + p) V - UT X y - VTyy + Qy _ 

( 2 . 20 ) 


2.4 Initial and Boundary Conditions 


For the final solution of the flow problem we need some initial values of the flow 
parameters to start with at time t = 0, the values are called initial conditions. Except 
as noted below, properties at each grid point are initialized as their respective free 
stream values. These initial conditions are 

u = Vqc cos a 
v = —Voc sin a 
P ~ Poc 

T = (2.21) 

After specifying these artificial initial conditions, suitable boundary conditions are 
specified. The boundary conditions applied are as follows: 

1. At the leading edge of the plate, no-slip is enforced (i.e U( 0 ,o) = ^(o,o) = 0.0 and 
the temperature (T( 0 , o)) and pressure (p(o,o)) are assumed to be their respective 
free stream values. 

2. At the inflow boundary (except the leading edge) and upper boundaries of the 
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Figure 2.1: Application of Boundary Condition 

domain, the x-component of velocity, u, the y-component of velocity, v, temper- 
ature and pressure are assumed to be their free-stream values. 

3. At the surface of the plate, no slip is specified on velocity (u = v = 0.0). Temper- 
ature (except the leading edge) is assumed to equal to the wall temperature T w . 
Pressure, (except the leading edge) is calculated at the wall by extrapolating 
from the values at the two points (i = 1 and i = 2) above the surface, as follows: 

P(oj) = ( 2 P(ij)) ~P(2j) (2.22) 

All properties on the outflow boundary of the domain (not including i = 0 and 
i - I MAX ) are calculated on the basis of an extrapolation from the two interior 
points, at the same j location. It is done as follows: 

u {i,J MAX) = (2U(,- ><wa A--1)) - x_ 2) (2.23) 
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2.5 Calculation of Skin Friction Coefficient and Rate 


of Heat Transfer 


2.5.1 Skin Friction Coefficient 


From the value of u and v on the surface and at grid point just above the surface, 
shear stress r xy is calculated using relation 


fdu s 

Txy ~ r Uv, 

Now skin friction coefficient is calculated from 


c / 


•xy 




For comparison with the results of boundary-layer equations C/ is calculated form 


c f 


i xy 


\peU 2 e 


where p e and u e are density and velocity u at the boundary-layer, i.e. at location when 
velocity profile begins to be vertical. 


2.5.2 Rate of Heat Transfer 


From the values of T and k on the surface and at the grid points just above the surface, 
rate of heat transfer is calculated form relation 


Qw 


-k 


dT 
dyj\ 


Stanton number for comparison with the results of boundary-layer equations is calcu- 
lated bv the use of relation 


where r = \/Pr. 


c H 


<7u' 


p e U e Cp + T 2j 


2 

2r e 
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Chapter 3 


Numerical Formulation and 


Solution Procedure 


The finite difference method has been used in this work to solve the system of equa- 
tions presented in Chapter 2, namely, Equation (2.2), (2.12) and (2.17). The principle 
used in this method is the conversion of partial differential equations to difference 
equations by the discretization of governing equations on a finite grid. The values of 
the flow variables are evaluated at each node and after each line step. The choice of 
the grid is very important for the accuracy of the grid. 


3.1 Discretization of Derivatives 


Each term in governing partial differential Equations (2.2), (2.12) and (2.17) can be 
written in a difference form using tavlor series expansion approximation of variable 
around the ith node(as described in Appendix A), x axis is on the surface of the wedge 
and y axis on the inflow boundary. 
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Figure 3.1: Finite Difference Grid Representation 
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du. 

U t - Ui 

-i 

W'~ 

Vi - y-, 

-i 

du. 

«i+i - 

Ui- 1 

ay i = 

y.+i - 

y t -i 


(F orward dif ference) 

( Backward dif ference) 
( Central Difference) 

( Forward dif ference) 

(■ Backward difference ) 
( Central dif ference) 


22 




3.2 Grid Generation 


Since with the increase in the Mach number, shock layer thickness decreases and heat 
transfer, pressure distribution on the surface increases, it needs finer grids throughout 
the domain (if we use uniform grid) for higher Mach numbers. But large number of 
grid points in the domain increases the computational time. To overcome this problem, 
we can use non-uniform grid with small vertical spacing near the plate and relatively 
large vertical spacing away from the plate. 


Both uniform and non-uniform grids have been used to obtain the results. For 
non-uniform grid non-uniformity in grid intervals is employed in y direction only while 
uniform intervals are used in x direction. Fig(3.2) shows the uniform and non-uniform 
grids respectively. 


3.3 MacCormack Scheme 


Writing the governing equations in the form 


dU 

dt 


dE 

dx 


9F 

dy 


(3.1) 


By means of Taylor series expansion the flow field variables are advanced at each grid 
point (ij) in steps of time as 


U« a ‘=UV+(^) At (3.2) 

__ \ / at* 

where once again U is a flow field variable (from the governing equations) assumed 
known at time t, either from initial conditions or as a result from the previous iteration 
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Figure 3.2: Uniform and Non-Uniform Grid Geometry 
in time.(^r) is defined as 



To obtain a value of (^r) the following steps are followed. 


( 3 . 3 ) 


1. (tet). . is calculated using forward spatial differences on the right hand side of 
v / * d 


the governing equations from the known flow field at time t. 
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2. From Step 1, predicted values of the flow-field variables (denoted by bar) can be 
obtained at time t + At as follows: 

US a ‘ = U‘,+ (^A( (3.4) 

Combining Stepl and Step2, predicted values are determined as follows: 

= U b - £ Ku - K) - £ (FU - Fb) < 3 - 5 ) 

3. Using rearward spatial differences, the predicted values (from Step 2) are inserted 

f QTJ \ 

into the governing equations such that the predicted time derivative J _ 
can be obtained. 

4. Finally substitute (7^)^ (from Step 3) into equation to obtain corrected 
second-order accurate values of U at time t + At. Hence combining Step3 and 
Step4 we get: 

= 5 [u*« + u:y - £ (e‘«< - Egfj) + ^ (F‘y< - 

(3.6) 


Steps 1 to 4 are repeated until the flow-field variables approach a steady-state value; 
which is the desired steady- state solution. 

In order to maintain second-order accuracy, the x-derivative terms appearing in E 
are different in the opposite direction to that used for ||, while the y derivative terms 
are approximated with central differences. Similarly, the v-derivative terms appearing 
in F are differenced in opposite direction to that used for |^, while the x-derivative 
terms in F are approximated by central differences. Hence for predictor step terms 
appearing in E are discretized as 


T xx{ij) — 


U {ij~ 1) 
Ax 


2 

3 Hi ' j) 


U {jJ) % +lj) 

Ax 2A y 


T xy — P(iJ) 


2A y Ax 
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. T(ij) T(i t j-i) 

q ‘ ~ ~ k ™ Tx 


and those in F are discretized as follows: 


'yy(ij) 




y (»j) 

2 

U (iJ+l) ^(tj-1) | ^(ij) 

Ax 


2Ax Ay 


T xy ~ H{iJ) 


u (*j) , ^oy+2) 

Ay 2Ax 


Qy ~ 


T M ~ Tq-ij) 

Ay 


while for corrector step terms appearing in E as follows: 
Txx (ij) = 2^(1, j) 


U (ij + 1) W (*\j) 

2 

u {tj+l) u (iJ) u (t+lj) 

Ax 


Ax 2Ay 


T *y ~ ^(ij) 


^(2+l,j) U (»J + 1) ^(»j) 

2Ay Ax 


Qx ^ 


(»J)' 




Ax 


and those appearing in F are discretized as follows: 


'mi j) 


^(tj) 


l ’(i+lJ) ~ V (’J) 


Ay 


Mu) 


^(i-nj) r (tj) 

2Ax Ay 


T xy ~ 


u {j+ hi) }Huj) , ^(tp+i) 

Ay 2Ax 
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Py — 


Ay 


After each predictor or corrector step, the physical variables are obtained by decoding 
the vector U; 


Pi,j = u i 

U 2 

Ui ’> = UT 

= U 3 
u,j Uj 

Et,, = U 4 

- _ U 4 v-h + v h 

Ui 2 


(3.7) 


with and e :j obtained, the remaining flow-field properties are determined 

by using the following relations: 



PiJ ~ PijRTij 


Pi,j = Po 


{ Tjj \ 2 Ip + 110 

V Tq ) T{j + 110 


k 


ij ~ 


Pr 


(3.8) 


where Pr is the Prandtl number taken to be constant at 0.71. 


3.4 Discretization of Boundary Conditions 


The y derivatives on the inflow boundary are obtained by forward difference while 
that on outflow boundary is determined by backward difference. Oh the other hand 
the x derivatives on the plate surface is determined by forward difference and that 
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on the upper boundary is obtained by backward difference, e.g. for plate surface and 


upper boundary 



U (2J) ~ 

Ay 


(-) 

\dy J ( IMAXJ ) 

and for inflow and outflow boundary 


u {IMAXj) ~ U(IMAX-lj) 

Ay 


^(i,2) ^(t, X) 

V ax /(i,i) _ Ax 





}Hi£MAX) ~ VjjJAfAX) 

Ax 


3.5 Calculation of Time-Step 


For the calculation of size of the time-step, the Courant-Friedrichs-Lewy (CFL) crite- 
rion is used [40]. aij is the local speed of sound- in metres per second and Kf is the 
Courant number. K f acts as a "fudge factor”, to make the solution stable. 

1 


(AtcFL)ij = 
where 


u. 


1,3 \ 


Ax 


+ 


M 

Ay 


+ G, 




v = max 


1 1 

T — 5 - + “7 — ^ + 2v 

Ax 2 Ay 2 

r |M»j (iWj/Pr) 

Pi -3 


— -—l 1 

Ax 2 Ay 2 1 


-1 


Now time step 

for 0.5 < Kf < 0.8 


At = min [Kj ( At CF L)ij 


(3.9) 

(3.10) 

(3.11) 


3.6 Solution Procedure 


As initial estimates, values are assigned to initialize the variables. MacCormack 
scheme is then applied to find the values of the variables at the next time step. Bound- 
ary condition is applied at the boundary points. Values at the next time is obtained 
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until the difference in density at any point in present time and that at the same point 
in the next time step is of the order of 10 -8 kg/m z . 

Shear stress and heat transfer on the plate surface are obtained by forward differ- 
ence. The solution procedure can be summarized as follow r s: 

1. Initiialize the variables u,v,p,T at all the grid points. 

2. Find other variables p, E, e, p, k using present value of u, v, p and T. 

3. Calculate the time step by the CFL criterion. 

4. Find the values of the variables in the next time step by MacCormack scheme. 

5. Apply boundary conditions. 

6. Check the convergence, if convergence condition is met go to Step 7. Otherwise 
go back to Step 2. 

7. Calculate skin friction and heat transfer on the surface. 
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Chapter 4 


Results and Discussion 


The numerical procedure given in the preceding chapter is implemented successfully 
for the Mach numbers 4, 8 and 12 over a flat plate and wedge at various angles. The 
computer system used for this purpose is HP-9000 of IIT Kanpur computer centre. 
On a pentium-II processor with 64 MB RAM the program needs 55 KB RAM and 
at 1009c CPU the program converges in 25 minutes. The disk space available to run 
the program is 20 MB. The results obtained are for surface length of 0.00001 m. and 
other properties are equal to those of standard sea level. The results obtained in this 
thesis for the hypersonic viscous flow regime are validated using results obtained in 
Ref. [40] . 

4.1 Validation of the Code 


The validity of the code developed is shown by obtaining the results of Mach number 
4 for 70X70 grid and comparing them with the results obtained in Ref. [40]. These 
comparison plots are given in Fig. (4.1). The results obtained for comparison are 
pressure, temperature and velocity profile, Mach number at the outflow boundary 
and normailized pressure distribution on the surface of the plate. The normalized y 
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distance ( y n ) has been obtained from 


where Re x 





4.2 Results for Flat Plate 


The plots shown in Fig. (4.2), (4.3) and (4.4) are of normalized pressure, velocity and 
temperature profiles, surface pressure distribution, skin friction coefficient and rate of 
heat transfer for various number of grid points (20X20, 50X50, 70X70 and 90X90) in 
the computational domain at Mach numbers 4. 8 and 12 respectively. Plots (A), (B) 
and (C) in all the three figures show that pressure, temperature and velocity profiles 
are approximately the same for grid 70X70 and 90X90, but plots (D), (E) and (F) 
show that normalized pressure, skin friction coefficient and rate of heat transfer on 
the plate surface at grid 70X70 and 90X90 are quite different. This observation shows 
that more accurate results can be obtained by taking finer grids near the surface. 

Fig. (4.5) shows that for Mach number 8 results, from 90X90 uniform and 40X40 
non-uniform grids are nearly the same, while they are very much different for Mach 4 
and 12. This show's that non-uniformity distribution of grids should be a function of 
Mach number. 

Plots showm in Fig. (4.6) are of normalized pressure and temperature profiles 
obtained at various times for Mach numbers 4, 8 and 12. These plots reveal that in 
shock capturing technique shock develops at any location wfith time instead of moving 
from one location to other. Pressure, temperature and velocity profiles obtained in 
plots A, B and C of Fig. (4.7) show reduction in shock layer thickness and increase 
in strength of the shock with increase in Mach number. Plots D and F show increase 
in pressure and heat transfer w'hile plot E show’s reduction in skin friction coefficient 
with Mach numbers. This reduction of skin friction coefficient does not mean that 
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frictional force reduces with Mach number, actually the (V^,) term in relation 


reduces C/ with Mach number. 



( 4 . 1 ) 


From Fig. (4.8) it is observed that Cfy/Re ~ x and c# \/Re x calculated from full 
Navier-Stokes equations are less than those calculated from boundary-layer equations. 
The reason for this is that the density of the flow at the boundary layer increases due 
to shock. Also it goes on increasing rapidly with Mach number. It is clear from these 
expressions 



T X y I p e U € X 

\PeV-l V Me 


(4.2) 


Cj/1 




(4.3) 


= 7-*=-, 

peUeCj, [T e -f r - T u .) V Me 
that the term ^fpl at the denominator increases with Mach number and hence de- 
creases the value of the term on the left hand srde of the Equations (4.2) and (4,3). 


4.3 Results for Wedges 


Plot (A) in Fig (4.9) shows coefficient of pressure distribution on a wedge surface 
for various Mach numbers. It is evident from the plot that as the Mach number 
increases, the coefficient of pressure decreases and comes closer to Newtonian pressure 
distribution. This reduction of coefficient of pressure does not mean that the actual 
pressure on the wall decreases with Mach number A/j n /fy • It is clear from the expression 


V Poc __ P Foe 

^ - ~ bpoo 


(4.4) 


that although Cj, decreases gradually with Mach number, pressure on the surface in- 
creases considerably as Mach number increases due to the variation. Plot (B) 

shows increase in pressure coefficient with wedge angle. It is also evident from plot 
(B) that deviation from Newtonian pressure distribution increases with Wedge angle. 
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Fig. (4.10) shows reduction in skin friction coefficient and increase in rate of heat 
transfer with Mach number for a wedge of wedge angle 10 degrees. Fig. (4.11) shows 
reduction in shock layer thickness and increase in it’s strength with Mach number 
for a Wedge at angle 10 degrees. Fig. (4.12), (4.13) and (4.14) show comparision of 
pressure, temperature and velocity profiles at the outflow boundary for a wedge at 
angle 10 degrees with those for a flat plate at Mach numbers 4, 8 and 12 respectively. 
The comperision of these plots with one another show that the rate of increase in 
shock strength and reduction in shock layer thickness with Mach number for a wedge 
is higher than those for a flat plate. 
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(A) Normalized Pressure Profile (C) Normalized Velocity Profile (E) Mach Number Profile 

(B) Normalized Temperature Profile (D) Normalized Surface pressure Distribution 

Figure 4.1: Comparison Between Computed Results and Results Given in Ref. [40] 
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Figure 4.4: Comparison of Various Parameters for Different Number of Grid Points 
at Mach Number 12 and Sea-Level Condition 
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(C) Skin Friction Coefficient at Mach 8 
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Figure 4.5: Comparison of Skin Friction and Rate of Heat Transfer on the Plate 
Surface for Uniform and Non-Uniform Grids at Different Mach Numbers 
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Figure 4.7: Comparison of Various Parameters for Different Mach Numbers 
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Figure 4.8: Comparison of Skin Friction Coefficients and Stanton Numbers for Bound- 
ary Lavers and Computed Results 
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(A) Comparision of Coefficient of Pressure for Various Mach Numbers 

(B) Comparision of Coefficient of Pressure for Various Wedge Angles 
(straight lines are Newtonian approximation at given angle) 

Figure 4.9: Comparison Plots of Wedge 
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Figure 4.10: Comparison of Coefficient of Friction and Rate of Heat Transfer for 
Various Mach Numbers for a Wedge at Grid 70x70 
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Figure 4.11: Comparison of Pressure, Temperature and Velocity Profiles at Various 
Mach Numbers for Wedge of Angle 10 Degrees and Grid 70x70 
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Figure 4.12: Comparison of Pressure. Temperature and Velocity Profiles for Flat Plat 
and Wedge of Angle 10 Degrees at Mach Number 4 and Grid 70x70 
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Figure 4.13: Comparison of Pressure, Temperature and Velocity Profiles for Flat Plate 
and Wedge of Angle 10 Degrees at Mach Number 8 and Grid 70x70 












Chapter 5 


Conclusions and Suggestions for 
Future Work 


5.1 Conclusions 


In this work the hypersonic viscous flow over a flat plate and wedges has been studied 
numerically, by the solution of Navier-Stokes equations using time marching finite 
difference method and shock capturing technique. These in turn have been further 
processed numerically to estimate the values of the individual skin friction coefficient 
and rates of heat transfer for various Mach numbers. 

Results show that increase in number of grid point in the domain increases the 
accuracy of the results. As expected, results of same and even higher accuracy have 
been obtained with less number of grid points when grids are nonuniform with higher 
density nearer to the surface of the plate. The distribution of nonuniformity for certain 
order of accuracy is found to be a function of Mach number. 

Increase in intensity of shock layer and decrease in thickness of the shock layer with 
Mach number is obtained as per expectation. Results after various time steps show 
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the development of shock and gradual increase in coefficient and rate of heat transfer 
with time. It has been shown that skin friction coefficient and Stanton number for 
full Navier- Stokes equations are significantly less than those obtained by solving the 
boundaiy-laver equations. 

Results obtained for wedge show the increase in skin friction coefficient and rate 
of heat transfer with increase in wedge angle. 


5.2 Future Scope of the Work 


Considerable scope exists for extending this work. The code can be used for the flow 
over a blunt-body by using a body-fitted grid and transforming it into a rectangular 
grid. It can also be made useful for an axisymmetric and three dimensional flow. For 
a three dimensional flow' problem, the third spatial parameters of the Navier-Stokes 
equations are also considered. 

The problem finds its important application in the analysis of effect of aerother- 
modynamics of the flow r on the structural dynamics of the body. As observed at higher 
Mach numbers temperature and pressure on the body are very high. Since high tem- 
perature reduces the stiffness of the struct ure, both temperature and pressure at hyper- 
sonic speed causes large deflection on the body (wings of the aircraft). This interaction 
between aerothermodynamics and structural dynamics (called aerothermoelasticity) 
can cause strctural failure. Thus limit of the speed can exactly be defined by consid- 
ering this effect. Chemical non-equilibrium effects can also be included in the code to 
make it accurate in the high Mach number range. 


Appendix A 


A.l Forward Difference Scheme 


Consider the Taylor series expansion of u j+1 around the point i. 

Ui+1 = u { + (*.'+i "**) , d 7 u (x i+1 - Xif 

dx u II + w li 2! + '“ (AJ) 

If the distance between x i+1 anf j x . j s small, its square terms will be even smaller so 
that the second deri\ati\e and higher order terms can be neglected. 


On . 

a? 


W.41 - 


(A.2) 


*i4l 

This L called forward difference of first derivative. In this approximation the order of 
accuracy is first. 


For a regular grid, the distance between two consecutive points wall be constant. 


i.e. x t+ ] 


x i ~ x 


»-t 


Ax. Therefore the forward difference scheme of first 


derivative becomes 



A. 2 Backward Difference Scheme 


Similarly, the forward difference involves expansion of ttj_i around i. 


Ui - 1 =Ui~ 


du (xj — Xi~i) 
dx li 1! 


d 2 u 


(.r, - x t _]) 2 
2 ! ^ ' 


(A.4) 


Again, if the distance between Xj_i and Xi is small, its square will be very small. So 
neglecting second order and higher order derivatives, 


chi | _ Ui — Mj_i 
dx Xi “ X^—i 


(A.5) 


This is the so-called backward difference of first derivative which is accurate to within 
first order. 


dv_\ _ Ui ~ »,-i 
dx ,<_ Ax 


(A.6) 


A. 3 Central Difference Scheme 


Subtracting (A.5) from (A.l) we will get 


du (x : -x : _i) n dhi (^-x^j ) 3 

Ul+1 - a,.! = 2-| t + 2—1, 


(A. 7) 


Neglect second order and higher order derivative terms, we can get 


du _ j _ n 1+ i - Uj-i 

Ox" Xi+i-Xi-i 


(A.8) 


This is called the central difference of first derivative. Here the truncation of term is 
from third only. Therefore central difference approximation has accuracy of second 
order. For a regular grid, it will become 


On. _ u, 4 i - Hi- 1 

ax 1, “ 


2 Ax Hrlffp*.! I 


(A. 9} 
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